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Atomistic models of epoxy polymers were built in order to assess the effect of structure at 
the nanometer scale on the resulting bulk properties such as elastic modulus and thermal 
conductivity. Atomistic models of both bulk polymer and carbon nanotube polymer 
composites were built. For the bulk models, the effect of moisture content and temperature 
on the resulting elastic constants was calculated. A relatively consistent decrease in modulus 
was seen with increasing temperature. The dependence of modulus on moisture content was 
less consistent. This behavior was seen for two different epoxy systems, one containing a 
difunctional epoxy molecule and the other a tetrafunctional epoxy molecule. Both epoxy 
structures were crosslinked with diamine curing agents. Multifunctional properties were 
calculated with the nanocomposite models. Molecular dynamics simulation was used to 
estimate the interfacial thermal (Kapitza) resistance between the carbon nanotube and the 
surrounding epoxy matrix. These estimated values were used in a multiscale model in order 
to predict the thermal conductivity of a nanocomposite as a function of the nanometer scaled 
molecular structure. 


I. Introduction 

T here has been a great deal of research effort in recent years on the subject of nanostructured materials . 1 

Structure across a range of length scales has been of particular interest. Among materials under study are 
nanocomposites, in which a particle with at least one dimension on the order of a nanometer is embedded in a 
polymer matrix. Polymers themselves are ‘nanostructured’ since they exhibit relevant structural features down at 
the nanometer length scale. Polymers are high molecular weight macromolecules, which leads to a variety of 
structural features in the condensed phase. In addition, the presence of contaminants such as moisture will 
contribute to the molecular structure of the polymeric system. In this paper, results are presented on atomistic 
simulations of epoxy polymers which are highly crosslinked macromolecular structures. These simulations are used 
to calculate macroscopic properties such as mechanical elastic constants and thermal conductivity. This is done with 
the goal of gaining insight into the way in which nanometer scaled structure influences macroscopic properties. 
Atomistic simulations of polymers in which the elastic constants are calculated have been performed for some 
time . 2 ' 1 However, atomistic simulations of crosslinked polymers have only been developed in recent years . 4 There 
have not been many fully atomistic simulations of crosslinked polymers in which mechanical properties have been 
calculated as a function of molecular structure. Multiscale simulations of thermal properties of polymer 
nanocomposites have been studied previously . 5 ' 7 In the thermal part of this work, the focus is on the effect of 
functionalization of the carbon nanotube within an epoxy matrix. 
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II. Method 


A. Building atomistic structures 

Epoxies are crosslinked polymer materials which are formed from the reaction of molecules which have epoxide 
groups and crosslinker molecules. In this paper, two types of epoxide containing molecules are used in simulations, 
a difunctional epoxy and a tetrafunctional epoxy. Figure 1 shows the two different epoxy types. In Figure la, the 
molecular structure of the diglycidyl ether of bisphenol A (DGEBA) is depicted. In Figure lb, the structure of 
tetraglycidyl diamino diphenyl methane (TGDDM) is depicted. The tetrafunctional epoxy contains four epoxy 
groups, leading to a potentially higher number of crosslinks per unit volume than the difunctional epoxy. In Figure 
2, the molecular structures of the crosslinkers used in the simulation are depicted. In Figures 2a-b, the structures of 
the two isomers of diethylene tolune diamine (DETDA) are depicted. In Figure 2c, the structure of diamino 
diphenyl sulfone (DDS) is depicted. Each of these crosslinkers contains two amine groups which can react with the 
epoxide. In some simulations, water (H 2 0) is also included in the simulation as a contaminant. In addition, a 
carbon nanotube (CNT) is also simulated with an epoxy in some of the simulations. The nanotube is functionalized 
by the chemical bonding of a functional group to the surface of the CNT. This functional group is an amide which 
contains an amine group that can participate in the crosslinking reaction. Figure 3a shows the chemical reaction of a 
crosslinking molecule containing an amine group reacting with a molecule containing an epoxide functional group. 
Figure 3b illustrates the product of the reaction from Figure 3a reacting with a second epoxide group. These 
reactions lead to the formation of highly crosslinked macromolecular structures. A mixture of an epoxy and a 
crosslinker will react according to the reactions depicted in Figure 3 to produce the cured epoxy network. In 
simulations in which a CNT is present, the CNT is simulated as extending through the periodic boundary conditions 
(PBC) in one of the directions. Further details on the construction of these atomistic epoxy models and of these 
nanocomposites have been presented previously. 5 ' 8 In this paper, results are presented from calculations performed 
on three structural types. The first two contained only cured epoxy with varying amounts of water added. The first 
of these was a DGEBA-DETDA epoxy. The second of these was a TGDDM-DDS epoxy. The third type built was 
a DGEBA-DDS epoxy with a CNT. No water was added to the CNT-epoxy nanocomposite. 

B. Molecular Dynamics simulation details 

In simulating these atomistic structures, the LAMMPS 9 software was used to run these molecular mechanics 
(MM) and molecular dynamics (MD) simulations. The force field used was AMBER 10 ' 11 while the TIP3P 12 
parameters were used to simulate water. A bond decrement method used for the CVFF force field in the Materials 
Studio software was used to estimate partial charges. 13 The unreacted epoxy and crosslinker molecules (including 
the CNT if present) were simulated in the condensed phase prior to crosslinking. The epoxy and crosslinker are 
present in a stoichiometrically equivalent ratio. This condensed phase reactant mixture was simulated with MD for 
2xl0 5 time steps (200 ps) using a 1 fs time step. This duration of MD time step was used in all the MD simulations. 
A final equilibrated configuration was taken and used for the construction of the crosslinked epoxy. 

Using this final static configuration of atomic coordinates of the condensed phase reactant mixture, a crosslinked 
network is formed by the addition of chemical bonds consistent with the chemical reactions depicted in Figure 3. 
The epoxide groups react with the amine groups to form chemical bonds. This results in a complex network 
structure. The method used is essentially that employed by most researchers for constructing atomistic crosslinked 
structures. 4 ’ 14 This new crosslinked structure is then used for the MM/MD simulation. 

The newly created atomistic configuration must be carefully equilibrated since the initial configuration is 
initially far from equilibrium. A short initial energy minimization is applied to the structure using a conjugate 
gradient algorithm in the FAMMPS software with an energy delta tolerance of 10" 4 kcal/mol. This is followed by an 
isothermal-isobaric MD simulation at a temperature of 300 K and a pressure of 1 atm for 2xl0 5 time steps (200 ps). 
The network formation procedure can be repeated in order to increase the degree of crosslinking. The degree of 
crosslink, a is a fraction of the number of epoxide groups which have reacted to become part of the network. In 
this stoichiometrically equivalent mixture, the maximum value for a is 1, which corresponds to all of the epoxides 
reacting with an amine. Once the desired degree of crosslinking, a , is achieved, a more complete equilibration is 
applied. An energy minimization procedure is applied with a restriction of energy delta tolerance of 10" 4 kcal/mol. 
The structure was simulated with MD for 10 5 time steps (10 2 ps) starting at a high temperature under isothermal- 
isobaric conditions. The structures were then cooled in a stepwise manner to achieve equilibrated structures over a 
range of temperatures. In the simulation of epoxy with no CNT present, the structures had water introduced in 
small percentages. The water molecules were added randomly to the structures. Following the addition of water 
molecules, energy minimization was applied. The structures were then subsequently equilibrated with MD using a 
stepwise cooling procedure described in the previous paragraph. 
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In order to obtain additional configurations for running deformation simulations, the MD simulation was run for 
an additional 100 ps. The final trajectory at each temperature is saved for running deformation simulations. The 
MD simulation was repeated for another 100 ps in order to obtain 3 separate atomistic configurations with their 
associated velocities for a given degree of crosslinking (a), moisture content and temperature (7). Although 

an MD run of this length (10 5 time steps, 10 2 ps) does not result in drastically distinct configurations, the properties 
calculated from these configurations did show statistical variation. 


C. Calculation of mechanical properties 

The elastic constants are calculated by monitoring the stress tensor while deforming the periodic boundary 
conditions during a MD simulation. This is similar to an established method which used static models, 2 but here is 
applied in a constant volume-constant temperature (NVT) ensemble. The 6x6 elastic constant matrix C is 
determined by the partial derivatives of the stress tensor, T, with respect to the deformation, £, as indicated by 
Equation 1. 



Three uniaxial tension deformations and three shear deformations are applied to each model. An isotropic symmetry 
is assumed, and the constants are averaged over the various orientations. The deformations were applied in a 
continuous fashion at every time step at a constant rate with respect to the original dimensions. The periodic 
boundary cell was deformed in a positive direction by 5% over a simulation period of 10 2 ps (10 5 time steps). These 
are typical rates of deformations for MD simulations and are necessary due to the limitation of the MD simulation 
method. 15 ' 16 The stress tensor was recorded at 10 time step (10 fs) intervals leading to a collection of 10 4 data points 
which were then averaged in batches of 10 3 time step intervals (10 ps). A least squares fitting method is used to 
obtain the elastic constants, fi and A from Equations 2-6. 3 Here, the contracted notation for the C matrix is used. 
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The Young’s modulus, E , and shear modulus, G, can then calculated from Equations 7-8. 

r 3/1 + 2// 

E = u — 

A + fi 

G = /u 


(7) 

( 8 ) 


D. Estimating the interfacial thermal resistance and thermal conductivity 

A multiscale modeling approach is used to estimate the thermal conductivity of CNT nanocomposites. 
Atomistic simulation is used to predict the interfacial thermal resistance between the CNT and the surrounding 
epoxy matrix. The interfacial thermal resistance was calculated 6 ’ 7 ’ 17 by instantaneously heating the nanoparticle 
atoms to a temperature of 500 K and then monitoring the difference in temperature, AT, between the nanoparticle 
atoms and the matrix atoms during a constant energy simulation (NVE ensemble). These data are then plotted as the 
natural logarithm of AT vs. time yielding a slope which is the negative inverse of the characteristic decay time, r. 
Equation 9 then gives the interfacial thermal resistance, Rk , where c T /A T is the heat capacity per area of nanoparticle 
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surface. This ratio ( Cj/A t ) is calculated assuming a value of 0.71 J/gK for c T 1 x This results in a value of 5.6xl0' 4 
J/m z K for Cj/A t for the CNT. 


Rk 


T 

( c t /A t ) 


(9) 


These interfacial resistance values from the heat transfer simulations are then used in an analytical formula 19 to 
calculate the thermal conductivities of nanocomposites. While molecular modeling yields R K values, an analytical 
model is applied in order to calculate the thermal conductivity of the nanocomposite material. This approach is 
taken because of the range of length scales that are relevant in the properties that are calculated. The interfacial 
thermal resistance depends on nanometer length scaled features of the nanocomposite and as such it is more 
reasonable to use a nanometer scaled modeling approach (MD). Thermal conductivity is a macroscopic property 
and a model is required which links nanometer scaled information (interfacial thermal resistance) with large length 
scaled features. This model assumes a random orientation of the dispersed nanotubes within the polymer matrix and 
uses an effective medium approach. This approach yields a formula for calculating the thermal conductivities of 
CNT composites. 19 This is summarized by Equations 10-13. 
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Material properties K e , K c , and K,„ are the thermal conductivities of the nanocomposite, the CNT and the polymer 
matrix, respectively. In evaluating Equations 10-13, a value of 0.2 W/mK is used for K,„. A value of 6000 W/mK is 
used for K c . 20 The diameter of the (8,8) SWCNT, d, is 1.077 nm, L is the nanotube length (assumed to be 10' 4 m), 
while/is the volume fraction, and a K is the Kapitza radius. 

III. Results 

A. Mechanical properties of the epoxy simulations 

Figure 4 shows the Young’s modulus for the TGDDM-DDS epoxy as a function of temperature for three 
different water contents. Figure 5 shows the Young’s modulus for the DGEBA-DETDA epoxy as a function of 
temperature for three different water contents. In both cases, the modulus decreases consistently with increasing 
temperature. The modulus of the TGDDM-DDS epoxy is much higher than that of the DGEBA-DETDA at 
comparable values a, f water and T. However, since (X is defined in terms of the degree to which the crosslinking 
chemical reaction proceeds and the TGDDM molecule contains twice as many epoxide groups as the DGEBA 
molecule, this parameter (a) is a bit misleading. The TGDDM-DDS model (a = 0.875) has density of about 1.8 
crosslinks/nm 3 , whereas the DGEBA-DETDA ( a = 0.86) model has a density of about 1.3 crosslinks/nm 3 . The 
increased number of actual crosslinks per unit volume is likely an important factor in the higher modulus of the 
tetrafunctional cured epoxy (TDDDM-DDS) versus the difunctional cured epoxy (DGEBA-DETDA). Figure 6 
shows the Young’s modulus as a function of water fraction f wa , er (in wt. %) for the TGDDM-DDS epoxy model 
while Figure 7 shows the Young’s modulus as a function of f water for the DGEBA-DETDA system. Here, the 
behavior does not seem quite as consistent as far as the Young’s modulus having a simple dependence on the water 
content. Although in general the moduli tend to decrease with increasing water content at the higher temperatures, 
there seems to be cases where the modulus actually appears to increase with increased water content. We have to 
consider several possibilities. One is that under certain conditions water may act as a plasticizer resulting in a 
decrease in modulus, and yet under other conditions, it may act as an anti-plasticizer causing an increase in modulus. 
Both conditions have been cited from experimental measurements in the literature. 21 The possibility exists that 
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under the high strain rates inherent in the MD simulation process, the relaxation mechanisms are different, resulting 
in discrepancies. There have been detailed analyses of the water-polymer matrix structure offered which suggest 
roles that various features such as free volume and the role of water molecules within the matrix have on these 
properties. 22 However, the possibility exists as well that these results are simply a matter of insufficient statistical 
sampling in these models. 22 In any case, we do not attempt to analyze this situation further within the current text. 

The effect that the degree of crosslinking has on volumetric properties was examined for the TGDDM-DDS 
epoxy. Figure 8 shows the specific volume for a TGDDM-DDS epoxy model as a function of temperature with a = 
0.56. Figure 9 shows the specific volume for a TGDDM-DDS epoxy as a function of temperature with a= 0.875. 
The lower a value in Figure 8 clearly shows greater expansion of the specific volume with increasing temperature. 
This aspect is especially noticeable at higher temperatures while at lower temperatures, the specific volumes 
converge to similar levels. 

B. Thermal conductivity 

Since we are interested in multifunctional aspects of materials, another property of interest is thermal 
conductivity. The dependence of thermal conductivity on nanometer scaled structural details was investigated by a 
multiscale modeling approach. MD simulation was used to calculate the interfacial thermal resistance between a 
CNT and the surrounding epoxy matrix as a function of the grafting density of covalently bonded amide functional 
groups. These amide functional groups participate in the crosslinking chemical reactions shown in Figure 3 such 
that the CNT becomes part of the crosslinked structure. Figure 10 shows the interfacial thermal resistance as a 
function of grafting density for the CNT-epoxy composite. Results are similar to those attained previously for an 
alkane chain group functionalized to a CNT embedded in a linear polymer matrix. 5 The interfacial thermal 
resistance declines noticeably for a certain range of grafting density (-0.008 A' 2 ) and then appears to decline very 
little with further increase in the grafting density. This result has some importance, as current synthetic technology 
may be limited in the amount of grafting density which can be achieved. Furthermore, a high value of the grafting 
density may result in degradation of various properties of the CNT. Using Equations 10-13, values are predicted for 
thermal conductivity as a function of grafting density for the CNT for several volumes of CNT loading in the 
nanocomposite. These calculated thermal conductivity values are plotted in Figure 1 1 . Fairly large increases in the 
thermal conductivity are predicted for the nanocomposite with moderate loads (5 volume %) of CNT. The strong 
effect of the functionalization of the CNT surface is seen as the grafting density of this functionalization is 
increased. 


IV. Conclusions 

Structural details at the nanometer length scale can influence macroscopic properties. The results here indicate 
how atomistic simulation can yield insight into these structure-property relationships. The elastic modulus was 
calculated for several epoxy systems. The results indicate a dependence of Young’s modulus on epoxy matrix 
molecular structure (difunctional versus tetrafunctional epoxy monomers, degree of crosslinking), on temperature 
and on the amount of moisture present. The dependence on moisture content remains somewhat elusive due to the 
complex issues involved. Speculation on these issues in the past has included free volume, moisture-matrix 
interactions, and statistical issues associated with the MD simulation. 22 Atomistic simulation, combined with a 
multi-scale modeling approach makes strong predictions for the effect of functionalization on the thermal 
conductivity of CNT based nanocomposites. 
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Figure 1 . The two epoxy molecules used in the simulation, a) the diglycidyl ether of bisphenol A (DGEBA) b) 
tetraglycidyl diamino diphenyl methane (TGDDM) . 
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Figure 2. The crosslinker molecules used in the simulation. a)-b) the two isomers of 
(DETDA). c) diamino diphenyl sulfone (DDS). 
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Figure 3. The chemical reactions of an epoxy and an amine. The upper chemical reaction (a) shows the reaction of 
an amine in a crosslinker molecule reacting with an epoxide group. The lower reaction (b) shows the amine in the 
product of the upper reaction undergoing a second reaction with another epoxide group resulting in a crosslink. 
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Figure 4. The Young’s modulus as a function of temperature for the TGDDM-DDS epoxy with degree of crosslink, 
a= 0.875, with varying water content. 
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Figure 5. The Young’s modulus as a function of temperature for the DGEBA-DETDA epoxy with degree of 
crosslink, a= 0.86, with varying water content. 
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Figure 6. The Young’s modulus as a function of water fraction, f wate . r , for the TGDDM-DDS epoxy with degree of 
crosslink, a= 0.875. The curves shown are for 3 different temperatures, 7 = 300 K, 350 K and 400 K. 
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Figure 7. The Young’s modulus as a function of water fraction, f water , for the DGEBA-DETDA epoxy with degree 
of crosslink, a= 0.86. The curves shown are for 3 different temperatures, T = 300 K, 350 K and 400 K. 
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Figure 8. The specific volume as a function of temperature for the TGDDM-DDS epoxy with degree of crosslink, a 
= 0.56. The three curves are for models with varying water content, f water = 0, 4.33 and 8.31 wt. %. 
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Figure 9. The specific volume as a function of temperature for the TGDDM-DDS epoxy with degree of crosslink, a 
= 0.875 . The three curves are for models with varying water content, f later = 0, 4.33 and 8.31 wt. %. 
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Figure 10. The interfacial thermal resistance, R K , between a carbon nanotube and the surrounding matrix as a 
function of the grafting density of functional groups on the surface of the nanotube. 
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Figure 11. The thermal conductivity of a carbon nanotube-epoxy nanocomposite as a function of the grafting 
density of functional groups on the surface of the nanotube. The 5 different curves correspond to different values of 
the nanotube volume fraction assumed. 
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